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I. INTRODUCTION 


Many structural concepts for large spacecraft applications involve 
a tensioned membrane surface. Since much of the structural weight in such 
spacecraft is associated with compression members necessary to equilibrate 
the membrane tension, a high premium is placed on designing structures 
with extremely small membrane tension. As a result, the elastic strains 
in the membrane may well be much smaller than anticipated thermal strains, 
so that a high likelihood exists for the development of wrinkled and slack 
regions within an otherwise taut membrane surface. The existence and 
severity of such wrinkled regions may have an adverse effect on the 
overall elastic stability of the spacecraft, as well as on the intended 
spacecraft performance if the membrane surface acts as a reflector with 
stringent requirements for geometrical accuracy of the surface. Thus, 
the problem of predicting the stresses and displacements within a partly 
wrinkled membrane surface is one of some current technological interest 
in the aerospace industry. 

In spite of the importance of the mechanics of wrinkling behavior of 
membranes the field remains largely unexplored. Apparently the earliest 
investigation in the field was reported more than 50 years ago by Wagner 
[1] who conceived "tension field theory" in order to explain the behavior 
of thin metal webs in beams and spars carrying a shear load well in excess 
of the initial buckling value. Wagner's method of analysis was based on 
lengthy geometrical considerations. Reissner [2] and independently Kondo 
[3] developed a simpler analysis based on straight-forward calculus, and 
presented the first exact solutions to problems involving a non-repetitive 
pattern of tension rays. Ia i [4] developed an analysis procedure based on 
a principle of maximum strain energy under given boundary displacements. 



This, and subsequent Japanese work on tension fields is well documented in 
a review paper (in English) by Kondo, Iai, Moriguti and Murasaki [5]. 

In a series of more recent papers, Mansfield [6-9] developed an analysis 
procedure which combines the "tension ray" concepts of Reissner and Kondo 
and a principle of Maximum "tension" strain energy similar to that of Iai. 

The first of his papers [6] developed an analogy with inextensional plate 
deformation theory, and worked example problems involving shearing and 
lateral contraction of membrane strips, and torsion of an annular membrane. 
Later [7], this work was extended to consider problems of load transfer from 
an elastic rod bonded to a flat membrane strip. Experiments were also 
reported which confirm predictions of the theory. Next [8], the analysis 
was generalized to include anisotropic and nonlinear membrane behavior typical 
of woven and fiberous materials. It was shown that such nonlinearities tend 
to amplify stress concentrations at corners and at the ends of a cut. Finally 
[9], the curved winkle patterns within hanging slack membranes was analyzed 
for various shapes of membranes and various support conditions. The curved 
wrinkles were shown to be governed by a one-dimensional diffusion equation 
and an analogy with heat conduction in a slab was noted. 

All the published work just cited considers static stresses and small 

deformations within fully-wrinkled flat membranes. A generalization to 

arbritrarily large deformations, with particular concern for the behavior of 

stretching skin, was recently presented by Danielson and Natarajan [10]. 

More recently, Wu [11] considered membrane wrinkling in the neighborhood of 

a sutured hole. Then Wu and Canfield [12] presented a general finite plane- 
on alw sis 

stress ^for wrinkling of flat membranes. 



However, for applications in space structures, a particularly important 
class of problems is that of partly wrinkled membranes. Such membranes 
contain both wrinkled and taut regions. A general theory for partly wrinkled 
membranes was developed by Stein and Hedgepeth [13] some 20 years ago. Their 
approach is based on experimental observations which show that when wrinkles 
develop within a membrane parallel to, say, the x-direction, the associated 
overall contraction in the y-direction exceeds that predicted by the Poisson's 
ratio effect. The additional average normal strain in the y-direction may 
be regarded as an "average wrinkle strain". For purposes of simplified analysis, 
these geometric features of wrinkling were incorporated in Ref. 13 into a 
Hookean material model by appropriately increasing the local effective value 
of Poisson's ratio in wrinkled regions. This effective value of Poisson's 
ratio may be determined by imposing the approximation that the local state 
of stress in a wrinkled region is one of uniaxial tension. The resulting 
theory retains the simplicity of form of the linear governing equations of 
elasticity, with the additional feature that the material parameters are dependent 
upon the local state of strain. Comparisons between the predictions of this 
theory and experimental results for some simple configurations [13,14,15] 
show that very satisfactory results may be obtained. Furthermore, of the 
available theories for wrinkling of membranes, the Stein-Hedgepeth theory 
appears particularly promising for finite element implementation. Recently 
[16], the approach was used to construct an algorithm for finite element analysis 
of flat membranes which contain taut, wrinkled, and slack regions. 
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With regard to curved membranes containing wrinkled regions, apparently 
the earliest work was done by Taylor [17] more than 60 years ago, in an analysis 
of parachute shapes. Taylor analyzed the geometry and stresses within a para- 
chute formed from an initially flat circular membrane by imposing the condition 
of zero hoop stress, and ignoring the stretching of the membrane near the 
crown. An extension of Taylor's work was presented more recently for isoten- 
soid surfaces by Houtz [18] and Mikulas and Bohon [19]. 

An analysis of axisymmetric doubly curved membranes which are formed 
from an intially flat membrane was presented by Mikulas [20]. The analysis 
allows stretching of the membrane surface near the crown to remove wrinkles, 
and includes a wrinkled region near the outer edge. The taut region near 
the crown is analyzed by a nonlinear membrane theory [21] based on Sander's 
theory of thin shells [22] with the omission of the bending terms. The wrinkled 
region near the outer edge is analyzed by imposing the zero hoop stress condition 
introduced by Taylor. Application is made to three example problems, including 
the pressurization of a pleated flat circular membrane attached to a rigid 
circular rim, the stretching of an initially flat circular membrane over a 
doubly curved, axisymmetric rigid mandrel, and the very large deformation 
behavior of a pressurized membrane cylinder subjected to a radial line load. 

Although it is not the primary focus of the reported research, the work 
of Zak on wave propagation within and vibration of partly wrinkled membranes 
is noteworthy. In a series of papers [23-30] Zak considers the stress conditions 
for local buckling within a thin membrane, and develops a continuum model 
for the propagation of shock waves, and "snaps" within thin membrane surfaces. 

His continuum theory for in-plane motion ignores the kinetic erergy of the 
out-of-plane motion in wrinkled regions, and results in an in-plane equation 
of motion which displays variable mass characteristics. As a result, his 



theory predicts a "damping" out of in-plane vibrations as the 
kinetic energy of the in-plane motion is converted into out-of-plane 
motion through shocks. Examples include a proof that static wrinkle 
patterns in fully wrinkled flat sheets loaded only along their edges 
must be straight lines. Also, one dimensional examples involving shock 
wave propagation in a freely hanging string, and vibrations of a single- 
degree-of-freedom oscillator constrained by an inextensible string are 
presented. 

Presented in this report are the results of one year of research 
supported by the National Aeronautics and Space Administration under 
Grant Number NAG-1-235 regarding the finite element analysis of wrinkling 
membranes. The research was performed at the University of Southern 
California (USC) in the Department of Civil Engineering, with technical 
assistance from Dr. John M. Hedgepeth of Astro Research Corporation 
of Carpinteria, California. 

Chapter two of this report presents the results of an analysis 
and exact solution for a problem involving axisymmetric deformations 
of a shallow curved membrane. The problem was considered in order to 
provide an exact solution and "benchmark" for calibration of future 
numerical examples. 

Chapter three describes the implementation of the numerical algorithm 
recently developed by Miller and Hedgepeth [16] on the SAP VII finite 
element code at USC. Chapters four and five then present an evaluation 
of this algorithm. The evaluation is based on a comparison of analytical 
and numerical results for stresses and displacements in two benchmark 
problems involving a partly wrinkled flat membrane. These comparisons 
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reveal a high degree of accuracy for the finite element algorithm. Further- 
more, convergence of the required iterative procedure in this nonlinear 
problem was achieved without excessive computation. 

I 
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II. AXISYMMETRIC DEFORMATION OF SHALLOW MEMBRANE 


2.1 General Analysis 


4 


% 


r 


Consider a shallow pressurized membrane of spherical radius a. For 
axisymmetric loading, the equilibrium equations are 



where N r , N 0 , and N rQ are the usual stress resultants, r is the radial 
dimension, p is the internal pressure, and P is an added center vertical 
load. 

The strain-displacement relations are 


e 

r 



e 


6 



( 2 ) 


Y r0 - *fe(?) 
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where u and v are the radial and circumferential tangential displacements 
and w is the normal displacement. We have assumed small strains and 
small slopes in comparison to unity except for the second-order term 
in the expression for e r ; this is necessary to be consistent with the 
displacement-dependent term in the third of Eqs. (1). Also, Eqs. (1) 
and (2) are valid only for r«a inasmuch as they apply to "shallow" 
configurations. 

We are interested in solving the equations for an inextensional 
membrane both wrinkled and unwrinkled. In an unwrinkled region, then 

£ r ' e 0 * ^re = 0 

In addition, in order that the principal stresses be non-negative, the 
stresses must obey the inequality 

N r N 0 » "re 

On the other hand, in a wrinkled region, 

N r N e= N re 
and 

fr _ N 0 

e e ' ^ 

The latter condition arises from the facts that the principal tension is 
to the wrinkles, and the principal strain is perpendicular to them. 


% 




parallel 
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Consider the annular region between r = and r = r^. Let w = u = v 
at r = r^, and let the load P and moment Q be applied to the "rigid" plug 
\ occupying the center. For small P and Q, there will be no wrinkles; 

for large P and Q, the entire annulus will be wrinkled. For intermediate 

* 

values, the innter part of the annulus will be wrinkled, and the outer 
portion will be unwrinkled. 

In all cases, the second of Eqs. (1) can be integrated to yield 



where 



is the shear stress at the inner boundary. 

2.2 Unwrinkled Region 

For the unwrinkled region, the displacements must represent rigid-body 
motion. Thus, 


u = v = w = 0 
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Solving the third of Eqs. (1) gives 


TT = 1 + P 

r 7 


( 6 ) 


The first of Eq. (1) then yields 


= 1 • T 

p 


( 7 ) 


where the various quantities have been nondimensional ized as follows: 


r,e 


— N 
pa ”r,0 


P = P/{nr 0 2 p) 

Q = Q/Ur^ap) 


P = r/r. 


In order that the membrane be unwrinkled 

fj fj > fj^ 

re Ye 

or 


p 4 > P 2 + Q 2 


(8) 


Clearly, if the right-hand side is less than unity, then there are no 
wrinkles. If it is greater than ( r 1 / r g) 4 « then the entire annulus is wrinkled. 



11 


2.3 Wrinkled Region 


Since 


N e' 


N re 

TTT 


we have 


Bp K) ■ 


(p% r ) 


(9) 


Integrating gives 


i 2 „ J- 

r p 4 




do) 


where N is the radial stress resultant at the center plug, 
r 0 


The third of Eq. (1) is nondimensionalized to give 


dw _ „ + ¥ 

dp pTL 


(u) 


where 


-ft) i 
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Substituting for N r from Eq. (10) and integrating gives 


2 k 

; - ; + <> - 1 JA 

-2 —2'' 

»: + q 
r 0 1 

)p 2 - ? 

W 0 + 2 

3 

("r + ° 2 ' 

\ 0 ) 

|3/2 




v 


X 



1 N r 

+ qV + 3N 2 P + (2 + 3P)Q 2 + — 

/ 0 I J-2 

J 3(N + 

\ r O 


—2 


X 


(1 + 3P)N 2 
r. 


(3 + 3P)Q 2 


] 


(12) 


where w Q is the vertical displacement of the central plug. 
The first two of Eqs.(2) in dimensionless form are 



where 



u 
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and 


E r.e = ( r n) Er ’ e 


As previously mentioned, the stress and strain relations must be 
related in a wrinkled region by 


E e N e * e r N r 


\ 


1 


Substituting and multiplying by an integrating factor gives 



This can be expressed as 




14 


Substituting from Eq. (11) gives 


d_ 

dp 



2N 


2 —2 

, . < P ± II- 

P N r 


(14) 


Going back to Eq. (13) gives 



r 


In order that one of the principal strains be zero, the shear strain 
must be related to the direct strains as follows: 


i 





4e e Q 
r 0 


Therefore, from the dimensionless form of the third of Eq. (2) 



or 


d_ 

dp 



22 fe 
p 3 K 


s 



V 



in which 


(15) 
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Define 


Then 


Integrating 


where 

v 

' F(p) = 


G(p) = 



(16) 


N r = ^ * / P 2 " M 

P 


/f®\ - p 3 

, (P 2 + P) 2 

\N J 2N/p 2 - M 

—2 2 — 

L N Z (p^ - M) J 


gives 

IT * "ZJ { ^ 2 [ F (P e ) " F <P>] - [ G <P e > “ G(P>][ < 17) 


I (P 2 . M) 3/2 + M(p 2 - M) 1 / 2 


1,2 - 5/2 . 2P + 3M . 2 -3/2 

- (P - M) + (p - M) 


+ (P + M) (P + 3M) (p 2 - M) 1/2 - M(P + M) 2 (p 2 - M)" 1/2 


(18) 
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Note that we have used the boundary condition that the circumferential 
strain must be zero at the outder edge of the wrinkled region, p = p g . The 
circumferential strain must also be equal to zero at the central boundary. 

Setting e Q (1) = 0 and solving for N gives 

i 

(19a) 

If the annulus is fully wrinkled then p fi = pj. Then N can be easily 
calculated as a function of M and P. Eq. (16) can be applied to determine 
Q and N p . 

If, on the other hand, the membrane is only partly wrinkled, then 
continuity of the radial stress resultant at p = p g requires that 

_ p 2 + ? 

(19b) 

Setting the two expressions for N equal to each other allows the determination 

of acceptable combinations of p , M, and F. Then TT and finally the 

accompanying values of ([ and TT can be determined. 

r 0 

For simplicity, the following analysis applies to pure force and y 

pure torque. 
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2.4 Pure Load 


For Q = 0, we set N = N ro , and M = 0. The expressions for the various 
quantities in the wrinkled region become 


N 


N 

e — 

r P 


2 3 — 

— — o-i o-i p 

w = w. + --- ■ ~ - M - f (p - 1) 

3N N 


( 20 ) 


2N p 


+ (2P - N 2 ) ^ + P 2 (p 


-] 


At the interface, e Q = 0, p = p , 


and 


N = p + — 
K e p 
K e 


( 21 ) 


Solivng for N, substituting into the last of Eq. (20), and solving 
for "P gives 


p - P„ 


/2pl + 4p2 + 6p e + 3 
5(2p e + 1) 


( 22 ) 


The deflection w must also be zero at p = p g . Therefore, 


(P e - 1)' 


B Po + P P + 1 Pp + 1 


P 
N 


3N 


(23) 
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Substituting for P from Eq. (21) gives 



2 2 Pe + 1 \ 

3 N ) 


(24) 

i 


Equations (22) through (24) enable the determination of the manner in 
which Wq and N and the extent of the wrinkled area change as P is increased 
above unity. When the entire surface becomes wrinkled (p £ = p), then the 
following results are obtained. 


p > 



15P 2 + 


10 


N = 


(p 2 + p x + l) 


P + 3 


/ 4 ^ 3 i 

(Pi + Pi + Pj 


+ 1 


5 (p + p + 1) 


(25) 




3P + p 2 + p x + 1 
3N 


P 1 + 1’ 


(P x - 1) 


(26) 


/ 


Equations (22) to (26) can also be used for negative values of 
P less than minus one by changing the sign of the right-hand side of 
Eq. (22). 
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2.5 Pure Torque 


For F equal zero, the expressions for F ang G in Eq. (18) become 
F(p) 


^ ~ M (p 2 + 2M) 


O 6 + 2M0 4 + 8M 2 d 2 - 16M 3 
G (p) = K “ 


5/p 2 - M 


(27) 


Then, Eq. (19a) becomes 


N 2 « | 


/ 6 — 4 —2 2 —3 

+ 2Mp + 8M* p - 16M* 3 
' e e e 


)/l - M - 


(1 + 2M + 8M - 16M )/p - M 

e 


5M 3 )/p 2 - M 


[p^ + Mp“ - 2M 2 )/l - M - (1 + M - 2M 2 )* / p 2 - M 


3 4 . -.2 2 

e ’ e 


(28) 


If P e is on the outer edge, the quantity N can be calculated as a function of M. 

Then Q and N can be determined from 
0 


Q = N A 


N = N A - M 


(29) 


For values of Q slightly greater than one, the region is only partly 
wrinkled. Then, using Eq. (19b) and solving for M gives 


5 - p; i 


A) 


(30) 
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Substituting from Eq. (28) yields the allowable combinations of p g and 

H for partial wrinkling presented in Table 1. Also given are the associated 

values of I) and TT found from Eq. (29) 
r 0 

The torsional motion of the hub is of interest. The angle of rotation 
is equal to 

0 V( v 

u = 

r o 


J 


Or, by using Eq. (15) 


17 



(31) 


where, as for the strains. 


\ 


a = 



n 


Substituting from Eq. (17) and judiciously using Eq. (28) and integrating 
yields 


ft » -2- 
3N 3 


-3N 2 


(^T7 - /i - m y 


a — 7 —2 

0 + 4Mp - 8M 

e 


/p 2 - M 


( 32 ) 


^ 1 + 4M - 8M 2 
x fc) 


/l - M 


The values for ft for the partly wrinkled case are given in Table 1. 



2.6 Numerical Results 


The load-displacement relations are shown in Figure 1 for pure load. 

The partly wrinkled membrane exhibits a softening behavior as the load is 
increased. When the wrinkling reaches the rim, the resulting fully wrinkled 
membrane stiffens with the application of increased load. 

The same behavior is exhibited for pure torsion as seen from Figure 2. 
It is interesting to note that, although deflections of the membrane begin 
when the loading parameter (P or Q) exceeds unity, significant movement 
starts occurring at a value of about two. This phenomenon is similar to 
that observed for cylinders in bending and is especially pronounced for the 
torqued dish. 



III. IMPLEMENTATION OF A FINITE ELEMENT ALGORITHM 
FOR A FLAT MEMBRANE 


3.1 Finite Element Algorithm 

Stresses and deformations in flat membranes may be described within 
the context of plane stress theory. For the class of problems under 
consideration, three regimes of structural behavior are possible. First, 
the membrane may behave in a fully taut manner, in which both principal 
stresses are positive. In general, this will occur whenever 

Ej > 0 and C£ > - ve^ 

where 

e l = 7 B e x + £ y’ + V< Sc - e y> 2 + fjy] 

e 2 ’ 7 [(c x t y ) - V(= x - + Y xy] 

In Eqs. (31) and (32), Ej and represent the local principal strains 

which are determined from the load-dependent strains e , e , and y . 

a y Ay 

In regions where Eqs. (33) hold, the stresses in the membrane may 
be determined from the well known plane stress elasticity relations. 
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where 


and 


where E and v are the modulus of elasticity and Poisson's ratio of the 
membrane material, respectively. 

In other regions the membrane may behave in a fully slack manner. 

In general, this will occur whenever 

e 2 - e l - °* (38) 

In such regions it is clear that the corresponding elastic principal 
stresses would both be compressive. However, since the membrane cannot 
support compressive stresses, the membrane is actually stress-free in 
such regions. Mathematically, this may be expressed as 


2' = ( °X* v V 
2 = ( V V Y xy )T > 


(36) 


E 


1 V 0 

v 1 0 

0 0 (l-v)/2 


(37) 


where 


a = D $ t 



(39) 


( 40 ) 
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Finally, it is possible for the membrane to develop wrinkles. In 
this case the stress state is regarded as uniaxial, with the tensile 
stress aligned along the direction of the wrinkles. Wrinkled behavior 
will occur whenever 

> 0 and Eg < “ v (41) 

In such cases, the material supports a tensile stress parallel to the 
priricipal direction associated with Cj, but is stress-free in a direction 
orthogonal to Cj. Mathematically, the stresses may be expressed as 

o = D w e (42) 

where 

2(l+P) 0 Q 

D w = | 0 2(1-P) Q (43) 

. Q Q 1. 


and where 



(cj - e 2 ) (ej - e 2 ) 


Further discussion of the description of stress-strain behavior within 
a wrinkled region, and the Stein-Hedgepeth constitutive model which 
forms the basis for the formulation presented herein, may be obtained 
from Ref. 16. 
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It should be pointed out that the choice of a D w matrix consistent 
with the Stein-Hedgepeth wrinkle model may not be unique. There may, 
in fact, be an infinite number of elasticity matrices capable of producing 
the correct stress-strain relations in a wrinkled element. However, 
the matrix of Eqn. (43) provides a bounded and symmetric matrix which 
has been successfully implemented and tested as reported herein and 
in Ref. 16. 

The stress-strain matrices presented above are strain-dependent 
and must be updated after each load increment. However, all other aspects 
of the problem formulation are identical to the approach used to solve 
any nonlinear stress/deformation problem where load increments and 
iteration for equilibrium are required. Consequently, the algorithm 
described above may be installed with relatively little effort in a 
variety of general purpose finite element computer programs which are 
capable of nonlinear analysis. 

3.2 Computer Implementation on the SAP 7 Code at U5C 

The algorithm presented in Section 3.1 was installed in the SAP 7 
computer code at USC where it now resides among the library of available 
material behavior for plane stress analysis. 

A description of the theoretical approach based on the principle 
of virtual work which is used to formulate the governing equations in 
the SAP 7 code, and the programmed methods for solving the resulting 
nonlinear equations is presented in Appendix A. 
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Most of the required modifi cations in the rather large SAP 7 code 
occurred in the subroutine ELPAL. A listing of the revised ELPAL routine 
which incorporates the wrinkle algorithm is presented in Appendix B. 

A listing of the input data used to generate the numerical results 
for verification examples 1 and 2 (discussed later) are presented in 
Appendices C and D. 


-r 
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IV. PURE BENDING OF A STRETCHED RECtANGULAR MEMBRANE 
(FINITE ELEMENT VERIFICATION EXAMPLE NO. 1) 

4.1 Problem Description 

As a first example of a partly-wrinkled flat membrane, consider 
a rectangular membrane which is uniformly pretensioned with normal stress 
o Q in the y-direction and with axial load P = o Q th in the x-direction, 
as shown in Fig. 3. Note that h is the length of the sides subjected 
to force P, and t is the thickness of the membrane. After pretensioning, 
an in-plane bending moment M is applied along the edges shown. As M 
is increased, eventually a band of vertical wrinkles of length b forms 
along the lower edge of the membrane as the normal strain e in this 
region becomes compressive. 

The solution for the stress and displacement fields within the 
resulting partly-wrinkled membrane is not trival. For example, the 
extent of the wrinkled region is not related in any simple way to the 
extent of the compression region within a similarly loaded flat plate. 

4.2 Analytical Solution 

A complete analysis of the problem just described is presented 
in Ref. [13]. It is shown in this reference that the extent of the 
wrinkled region may be determined by 


b _ 3M 1 

TT " "PK ' 7 . 


( 45 ) 
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Furthermore, the overall moment-curvature relation for the membrane 
is given by 


(46) 

and the stress field within the membrane is given by 






* 

(49) 


These analytical results are used to evaluate the accuracy of the numerical 
solution generated by the finite element model discussed in the next, 
paragraph. 
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4,3 Finite Element Modeling Assumptions 

A finite element model of the problem described above was created 

‘ using the rectangular grid shown in Fig. 4. The model consists of fifty 

isoparametric quadrilateral elements, each of which contains four internal 

* 

integration points at which stresses are determined. Displacements 
at the four corners and midpoints of each side of each element are also 
determined. The grid contains a total of 181 of such nodal points. 

The number of unconstrained degrees of freedom in this model is 362. 

The problem shown in Fig. 3 is essentially one of specified loading 
conditions around the perimeter of the membrane. The specification 
is not complete, however, since only the resultants P and M of the stress 
distribution along the left and right edges of the membrane are specified, 
and not the detailed stress distribution itself. This situation is 
also unsatisfactory from the numerical modeling viewpoint because the 
model is not restrained from rigid body motion and, the resulting global 
stiffness matrix would be singular. Therefore, some external constraints 
on the displacement of at least two nodes is necessary in order to develop 
a satisfactory finite elment model. It is important to note that such 

* constraints will generally lead to local deviations in the stress and 
displacement fields from those predicted by the analytical solution. 

* 

After investigating several edge constraints and loading conditions, 
a model was adopted which consists of the constraints u x = 0 at each 
of the eleven nodes along the left edge of the membrane, where u is 

A 

the nodal displacement in the x-direction. All nodal points along this 
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edge are free to move in the y-di recti on except the node at the center 
of this edge, which is constrained such that u = 0. 

Along the upper and lower edges of the membrane, the loading conditions 
used in the finite element model consist of vertical tensile loads 
applied at each node. The magnitude of these nodal forces is determined 
such that they are equivalent to a uniform tensile normal stress of o Q 
along each edge. 

In order to apply the resultant tension force P in the x-direction 
and bending moment M, the nodes along the right edge of the membrane 
model were attached to a finite element model of a very stiff beam upon 
which external forces were applied. The attachment between the membrane 
and beam models was accomplished by requiring continuity of displacements 
in the x-direction at each node. However, displacements in the y-direction 
for nodes belonging to the membrane were not required to be the same 
as those for nodes belonging to the beam, except at the node in the 
center of the edge. In addition, nodal forces were applied to the beam 
in such manner that the resultant force was P and the resultant bending 
moment was M. 

Numerical results were generated by first applying the pretensioning 
forces. These forces correspond to P along the right edge and along 
the upper and lower edges. After the equilibrium configuration of the 
pretensioned membrane was obtained, a bending moment M was applied in 
small increments to the stiff beam along the right edge of the membrane. 

An iterative solution for equilibrium displacements was generally required 
after each increment of bending moment. 



4.4 Comparison of Finite Element and Analytical Results 


The qualitative nature of the reults of the finite element simulation 
are shown in Figs. 5 and 6. Shown in Fig. 5 are the directions of the 
edge displacements along the right edge of the membrane. (The displacements 
are not to scale.) This displacement pattern reveals the downward translation 
and clockwise rigid rotation which are expected for the edge displacements 
of a cantilevered beam subjected to a pure clockwise bending moment. 

Shown in Fig. 6 is a plot of the directions of the principal stresses 
within every element for the case of an advanced loading state. In this 
state the applied bending moment is so large that about 75% of the membrane 
surface is wrinkled. The principal stresses are shown centered on each 
internal integration point as orthogonally directed line segments. The 
line segments define the orientation of the principal axes of stress. 

Wrinkled regions are indicated by a single line segment in a direction 
parallel to the wrinkles. The length of the arrows is not proportional 
to the magnitude of the principal stress, and in this sense the figure 
is not to scale. The results indicate that the principal axes of stress 
are all aligned along the x and y coordinate axes, as expected. 

A comparison of the overall moment-curvature behavior of the membrane, 
as determined by the finite element model and by Eq. (46) is shown in 
Fig. 7. Note that the numerical results, even for very large curvatures, 
are very accurate. In fact, the errors are so small as to be nearly 
imperceptable at this scale. 
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Shown in Fig. 8 is a comparison of analytical (Eq.(45))and numerical 
results for the wrinkle band width, (b/h). Since this band width is 
not directly available from the finite element code, it was necessary 
to estimate the location of the boundary between wrinkled and taut regions 
by linear extrapolation of numerical results for o x . This was accomplished 
by plotting numerical values for a x vs. y (as in Fig. 9), then fitting 
a curve through these points and extrapolating to a = 0 in order to identify 

A 

the corresponding value of y = b. Repeating this process for three different 
stress states corresponding to three different levels of bending moment, 
the data were obtained for Fig. 8. Again it is seen that the errors 
in the numerical results are less than two percent. 

Shown in Fig. 9 are curves for o vs. y for three different levels 

A 

of bending moment. The analytical results were generated from Eq. (47), 

and the numerical results were obtained directly from the SAP 7 output 

files. Although they are not plotted, the numerical and analytical 

results for t agreed on the value of zero in every element, and a 

y 

was observed both numerically and analytically to be o Q in every element. 

The numerical results shown in Fig. 9 correspond to the o stresses 

A 

at each internal integration point within the vertical strip of five 
elements located just left of the center of the membrane model shown 
in Fig. 4. Results for an element strip along the extreme left or right 
edges of the membrane vary slightly from the plotted results due to 
imperfect boundary conditions in the model, as previously discussed. 

The accuracy of the numerically determined values for o x is again 
found to be very high at every location except for the special case 
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when the integration point happened to be located very near the boundary 
y=b. Since no efforts were made to devise an advanced algorithm for 
such boundary points, the discrepancy in this special case is not surprising. 
However, the principal observation from Fig. 9 is that the finite element 
model is capable of providing very accurate results for the detailed 
stress distribution within the membrane, at almost every point. 


V 
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V. PURE ROTATION OF A HUB ATTACHED TO A FLAT STRETCHED 
MEMBRANE (FINITE ELEMENT VERIFICATION EXAMPLE NO. 2) 

5.1 Problem Description 

As a second example of a partly-wrinkled flat membrane, consider 
a rigid circular hub of radius a attached to a flat stretched membrane 
as shown in Fig. 10. The membrane forms an infinite flat sheet and 
extends indefinitely in all directions. Let the rigid hub be perfectly 
bonded to the membrane before pretensioning. After the hub is attached, 
let the membrane be subjected to a uniform edge tension at infinity, 
so that the stress state in the membrane far from the hub is isotropic 
with principal stress o Q . After pretensioning, the hub is subjected 
to a pure twisting moment M t , as shown. For sufficiently large M t the 
membrane stress state consists of an exterior taut region and an annular 
interior region which is wrinkled. The extent of the wrinkled region 
is measured by the wrinkle radius R. 

5.2 Analytical Solution 

A detailed analysis of a class of problems very similar to the 
one just described is contained in Ref. [14]. In particular, from 
Appendix B, of Ref. [14] in the limiting case of an infinite membrane 
(b the wrinkle radius R corresponding to a prescribed twisting 
moment M t is governed by 


1 

I 


+ 


1 

F " 



0 


(50) 
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1 


where 



B = tf 2 ^- 1 

IT 2 


( 51 ) 



(52) 


(53) 


In deriving Eqs. (50) through (53) it was assumed that the value of 
Poisson's ration for the membrane material is v = -j, and that the thickness 
of the membrane is t. 

Furthermore, after correction of a typographical error in Eq. (37) 
of Ref. [14], the induced angle of twist <J> of the rigid hub is found 
to be governed by 


-r - 3M 
" T 



(54) 


where ^ = 2G4>/a Q and G is the shear modulus of the membrane. 

The simultaneous solution of Eqs. (50), (51), and (52) requires 
an iterative numerical approach. However, it can be shown that the 
response is linear (i.e., no wrinkles occur) for M < V 3/2 . In this 
regime, one finds that M = and R is not meaningful. The onset of 
wrinkling occurs at M =>[i/2 at which IT = 1. For M the corres- 

ponding nonlinear solutions are presented in Table 2. 
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Furthermore, the principal stresses in the membrane may be shown 


to be 


>1 . 


«Vr) 


u iz\ 2 


£4 - (M/r) 


1 + (R/r ) 2 


1 < r < TT 


R < r 


(55) 


(Tt/r ) 2 


1 < r < R 


R < r 


(56) 


where 


r = r/a 


(57) 


and r is the radial distance from the center of the hub. 
Also shown in Table 2 are numerical values for 


16 T r0 (? = 4) 
° r (r = 4) 


M 

i+ 4^ 2 -f 


(58) 


The values for k are used later in determining an equivalent value of 
M in the finite element model for this problem. After the appropriate 
K has been identified, the analytical values for ft, 4 >, ( o j/o 0 ) and (° 2 /° 0 ) 
are compared with the corresponding results from the finite element model, 
which is described in the next paragraph. 
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5.3 Finite Element Modeling Assumptions 

The creation of a finite element model for this infinite membrane was 
accomplished by first imagining the problem as being the superposition 
of an interior problem and an exterior problem which are pieced together 
at a common interface. Let the interface be a circle of arbitrarily 
chosen radius r=4a in Fig. 10. Then the exterior problem (r>4a) contains 
no wrinkled region and may be solved analytically by simple two-dimensional 
elasticity theory. The interior region so created (r<4a) is a finite 
dimensional annulus which contains a wrinkled region. 

The appropriate interface conditions at r=4a are continuity of o r 
and T r0 , and of displacements u p and u Q . The approach employed consisted 
of applying prescribed o r and i r9 along r=4a for both problems. The 
hub in the interior problem was then regarded as fixed against displacement, 
as was the outer edge in the exterior problem. The angle of twist in the 
corresponding composite problem was obtained by adding the angle of 
rotation between the rim r=4a and fixed hub in the interior problem, and 
the angle of rotation between the rim r=4a and the fixed boundary at 
infinity in the exterior problem. 

The appropriate values of M t and o Q were obtained from the known 
values of o r and T r0 applied at r=4a. This was accomplished for given 
applied stresses o f and i r6 at r=4a (which are sufficiently small that 
R<4a) by first computing k from Eq. (58). For a given value of k, the 
corresponding value of M may be obtained from Table 2 or from Eqs. (50), 
(51), (52), and (58) by an iterative solution. o 0 may then be found 
from the* equilibrium relation 



( 59 ) 
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where x r0 is known. Once M is determined for a given loading case in the 
finite element model, M t is determined as well as o Q , and the nondimensional 
stresses and displacements may be compared with the analytical values previously 
discussed. 

The finite element model of the interior problem was created using 
the quasi -circular grid shown in Fig. 11. The model consists of 36 isopara- 
metric quadrilateral elements, each of which contains four internal integration 
points at which stresses are determined. Displacements at the four corners 
and midpoints of each side of each element are also determined. The grid 
contains a total of 132 of such nodal points. The number of unconstraining 
degrees of freedom in this model is 264. 

The nodes along the hub-membrane interface were considered fixed against 
displacement. That is, u r = u 0 = O along r=a. Along the outer edge r=4a, o r 
and T rQ are prescribed independently, but are applied such that t =0 initially 
while o is increased to a constant pretension value. After the pretension 
in o r has been achieved, then t is increased in steps. As a result wrinkles 
are eventually initiated, and the wrinkle radius R increases with increasing 

T r0‘ 

The principal stresses Oj and at each integration point within each 
element are determined by the SAP finite element program, and may be read 
directly from the output files. The same is true of the nodal displacements 
u r and u 0 . The angle of twist for this interior problem may be computed 
from 


♦ 


I 


u 0 (r = 4) 


4a 


(60) 
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where u Q (r=4) represents the average tangential displacements around the 
outer perimeter r=4a of the finite element model, as read directly from 
the computer printout. 

The angle of twist <t>^ for the corresponding exterior problem may be 
derived from two-dimensional linear elasticity theory. The appropriate 
boundary conditions for the exterior problem are (1) a r and T rQ * r: are prescribed 
at r=4a (2) o r , o e -»-a 0 and u Q -+0 as r* °°. The solution for <f>£ may be obtained 
as 




(61) 


The total angle of twist for the numerical model is then defined as 

* = *E + *1 (62) 


5.4 Comparison of Finite Element and Analytical Results 

The qualitative nature of the numerical results for the interior problem 
are shown in Fig. 12. The directions of nodal displacements (not to scale) 
are shown for a typical one-quarter sector of the annular interior region. 

The displacements are seen to be primarily clockwise rotation, with a inward 
radial component noticeable near the hub. 

Shown in Fig. 13 is a comparison of analytical (eq. (55)) and numerical 
values for the maximum principal stress (a /cr ) as a function of radial 
position (r/a) for four different load cases. The numerical values for 
Oj were obtained directly from the SAP 7 output files for a typical group 
of these elements which form a sector of the grid in this axisymmetric problem. 
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The accuracy of the numerically determined values of a ^ is seen to be very 
high in all cases. 

A comparison of analytical (Eq. (56)) and numerical values for the 
minimum principal stress (o 2 /o 0 ) as a function of radial position (r/a) 
is shown in Fig. 14, for four different load cases. Again the numerical 
values for a 2 were read directly from the SAP 7 output files for elements 
in a typical sector of the grid. The accuracy of these minimum principal 
stresses is seen to be considerably less than that of the maximum principal 
stresses shown in Fig. 13, but nevertheless the numerical values may still 
be quite adequate for many engineering purposes. 

Shown in Fig. 15 is a comparison of the analytical (Eq. (50)) and numerical 
values for the wrinkle radius (R/a) as a function of the applied torque K. 

Since R is not directly available from the SAP 7 output files, it was neces- 
sary to estimate R by curve fitting the results for (o 2 /o Q ) vs. (R/a), and 
then extrapolating to the case (°2/ a 0 ) ~ This process was repeated for 
each load case to obtain the data plotted as numerical values in Fig. 15. 

The accuracy of the numerical values so obtained is found to be adequate 
for many engineering purposes. 

A comparison of the analytical (Eq. (54)) and numerical values for 
the overall angle of twist i as a function of the total applied torque M 
is shown in Fig. 16. The numerical values for 4> were obtained from Eqs. (60), 
(61), and (62) as previously described. As shown in the figure, the finite 
element model produces results of acceptable accuracy for small torques, 
but the errors increase as the torque levels increase to a maximum of about 
five percent at M = 5. As expected, the errors indicate that the finite 
element model, which in this example involves only three elements in the 
radial direction, is too stiff in an overall sense. It is expected that a 
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finite element model with more elements would result in more flexible 
behavior, and a reduced error in Fig. 16. 
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TABLE 1. PARAMETER COMBINATION FOR PARTIAL 
WRINKLING WITH PURE TORQUE 

P e H A S \ 0 * 

1.0 0.7071 1.4142 1.0000 1.0000 0.0000 

1.2 .6320 1.6020 1.2375 0.9718 .0004 

1.4 .7548 1.7854 1.5516 .8841 .0057 

1.6 .8570 1.9617 1.8161 .7418 .0299 

1.8 .9284 2.1310 2.0533 .5703 .0968 

2.0 .9681 2.2972 2.2603 .4102 .2325 

2.5 .9959 2.7267 2.7211 .1738 .9413 

3.0 .9993 3.1818 3.1807 .0839 2.2019 

3.5 .9998 3.6522 3.6519 ‘ .0456 4.0771 

4.0 1.0000 4.1312 4.1311 .0270 6.6473 

4.5 1.0000 4.6154 4.6154 .0171 10.0394 

5.0 1.0000 5.1031 5.1031 .0113 14.3224 

5.5 1.0000 4.4932 5.5932 .0078 19.6060 

6.0 1.0000 6.0851 6.0851 .0055 24.9899 
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TABLE 2 

Analytical Solutions for Stresses and Displacements 
in the Hub Rotation Example of Figure 10 


~w 

0.8683 

0.8776 

0.8891 

0.9816 

1.0976 

1.3311 

1.5656 

1.7987 

2.0286 

2.2545 

2.4765 

2.6949 

2.9101 

3.1226 

3.5409 

3.9523 

4.3585 

4.7607 

5.1597 

5.5563 

5.9510 

6.3440 

6.7356 

7.1262 

8.2929 

9.0677 

10.2268 

10.9981 


K 

1 

1.001 

1.0050 

1.0100 

1.0500 

1.1000 

1.2000 

1.3000 

1.4000 

1.5000 
1.6000 
1.7000 
1.8000 
1.9000 
2.0000 
2.2000 

2.4000 
2.6000 
2.8000 

3.0000 
3.2000 

3.4000 

3.6000 
3.8000 

4.0000 

4.5000 

5.0000 

5.6000 

6.0000 


♦ 

\[V2 

0.8683 

0.8776 

0.8892 

0.9849 

1.1123 

1.4019 

1.7454 

2.1403 

2.5753 

3.0386 

3.5216 

4.0190 

4.5276 

5.0455 

6.1045 

7.1901 

8.2993 

9.4301 

10.5810 

11.7507 

12.9385 

14.1432 

15.3642 

16.6007 

20.3968 

22.9943 

26.9814 

29.6953 


3 

3.0040 

3.0201 

3.0403 

3.2090 

3.4387 

3.9786 

4.6529 

5.4774 

6.4467 

7.5455 

8.7593 

10.0773 

11.4924 

13.0000 

16.2792 

19.9002 

23.8547 

28.1382 

32.7482 

37.6829 

42.9412 

48.5223 

54.4256 

60.6508 

81.2545 

96.5952 

122.0112 

140.5578 


k 

0.8398 

0.8420 

0.8510 

0.8622 

0.9517 

1.0637 

1.2869 

1.5057 

1.7152 

1.9122 

2.0957 

2.2656 

2.4225 

2.5673 

2.7007 

2.9355 

3.1321 

3.2945 

3.4266 

3.5319 

3.6137 

3.6748 

3.7179 

3.7456 

3.7599 
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Fig. 1 - Applied Torque vs. Angle of Twist 
for Axisymmetric Shallow Membrane 
Problem 



LOAD 
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Fig. 2 - Applied Load vs. Vertical Displacement 

for Axi symmetric Shallow Membrane Problem 
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Fig. 3 - Flat Stretched Membrane Subjected 
to Pure Bending Moment 





Fig. 4 - Finite Llement Model for the 
Membrane in Fig. 3 
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Fig. 5 - Qualitative Plot of Edge Displacements 
from Numerical Solution for Pure Bending 
of Rectangular Sheet. (Displacements 
not to scale.) 
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Fig. 6 - Qualitative Plot of Principal Stresses 

from Numerical Solution for Pure Bending 
of Rectangular Sheet (Advanced Loading State.) 
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Fig. 7 - Moment-Curvature Relation for 
Pure Bending of Rectangular 
Sheet 
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Ph 

Fig. 8 - Wrinkle Height vs. Bending 
Moment for Pure Bending of 
a Rectangular Sheet 


x 

h 

Fig. 9 - Maximum Principal Stress (o x ) vs. 
Vertical Position (y ) for Pure 
Bending of a Rectangular Sheet 
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Fig. 10 - Flat Stretched Infinite Membrane 

Subjected to a Pure Twist M Through 
an Attached Rigid Hub. 









Fig. 11 - Element and Node Configuration Fig. 12 - Qualitate Nodal Displacements 

for Finite Element Model of Hub from Numerical Solution of Hub 

Rotation Interior Problem Rotation Interior Problem. 

(Displacements not to scale.) 





r_ 

a 


Fig. 13 - Maximum Principal Stress vs. Radial Position for 
Four Load Cases (Hub Rotation Interior Problem) 
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$ 


Fig. 16 - Nondimensional Angle of Twists vs. 

Applied Torque M for Hub Rotation 
Composite Problem 
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APPENDIX A 


FORMULATION AND SOLUTION OF 
NONLINEAR STRUCTURAL PROBLEMS 
USING SAP 7 


This section will discuss the formulation of, the governing equations 
starting from the principle of virtual displacements and methods for solving 
the resulting nonlinear equations. The discussion of the formulation of the 
equations follows references [1] and [2], 


FORMULATION 


The motion of a general body is shown in Fig. 1. The configuration is 
known at times 0 and t and the objective is to determine it at time t + At. 

In the following derivation a left superscript indicates the time when the 
quantity occurs and a left subscript indicates the configuration with respect 
to which the quantity is measured. In the case of derivatives, a left sub- 
script indicates the time of the coordinate with respect to which the quantity 
is differentiated. Thus, 



FIG. 1 MOTION OF BODY IN CARTESIAN COORDINATE 
SYSTEM (FROM REFERENCE (1) ) 
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iu 

E t 

t*tt u U 8( ,4A, x.) 


( 1 ) 


All tensors are referred to Cartesian reference frames. The principle of 
virtual displacements, written for the current configuration (time = * + At ) 
is 




where 


„ I* ^^ tk)hUk f dA) + 


/ °pC Lt 0 f k ) fdV) 
°v 


The quantities ( t4AtT *>) are the Cartesian components of the Cauchy 
(true) stress tensor at time t + &t , and ( f4A ©f*) are surface tractions 
and body force components at time r+Ar but measured with respect to the 
configuration at time = 0. The variation bu k is an infinitesimal variation 
in the current displacement component (r + Ar,^) . The summation convention 
for repeated indices is used here. The variation in true strain corres- 
ponding to the infinitesimal variation in the displacement field is 


>! 


6 4 


1 (r + Ar„ 




’“V 


(3) 


In dynamic analysis the body force components include inertial effects. 

Since the configuration at the current time r+Af is unknown, the 
principle of virtual displacements, Eq. (1), must be expressed in a form 
in which all variables are referred to a known state. Then the integration 
will be performed over known volumes and areas. If the static and 
kinematic variables are referred to the initial state (time = 0), the 
formulation is known as Total Lagrangian (TL). If they referred to the 
previous known state (time = t), the formulation is known as the Updated 
Lagrangian (UL). The remainder of this derivation will follow the total 
Lagrangian formulation. References [1] and [2] present details of both the 
TL and UL formulations. 

In the TL formulation the principle of virtual displacements becomes: 




'*»'/> ( 4 ) 
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t*tt s 

in which •' a are the components of the 2nd (symmetric) Piola-Kirchoff 
stress tensor and are the components of the Green-Lagranqe strain 
tensor which is written in terms of the current displacements t * tt u k : 




+ t+Af u 
+ o u k,i 




(5) 


If incremental static and kinematic quantities are defined as 


o^= ,+ A, 5 ( y-'S ( y 


€.,• = 


= t * Af, 


o'- n 

U: St* Lt u 


ij 0*IJ 


> i 


( 6 ) 


the total Green-Lagrange strain increment can be decomposed into linear and 
nonlinear parts 


0 € ij " 0 e ij + »%i 


(7) 


with 


•*ii “ ^ \e u U*o u i,l * Q> u k.i) Wj> 

* (>*./» < 8 > 


In addition, the increments of 2nd Piola-Kirchoff stress tensor components 
can be related to the increments of Green-Lagrange strain tensor components 
by the linear constitutive law 


which is an approximation, since • c iin changes along the path in the 
finite increment At. Furthermore, the variation in the total strain 
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4 *•«/>* is equal to the variation in the total strain increment Hpty) • 
This equivalence and the relations (6), (7), and (9) can be used to express 
Eq. (4) in the form 



ijrt 



. t+*t R 


-/ O S iJ B (o e lj)° dV 
°V 


( 10 ) 


This variational principle represents a nonlinear equation for the incre- 
mental displacements u/ . The Updated Lagrangian formulation leads to a 
completely analogous expression in which the subscripts and superscripts 
(o) are replaced by (t) and the equivalence of the 2nd Piola-Kirchoff 
stress tensor t S u to the Caunchy (true) stress tensor t r t j is noted. Of 
course, the constitutive tensor for the U.L. formulation is now t c an , 
which takes on different values from 0 c Hn . The relationship between the 
two constitutive tensors is given by 


t C mnpq 1 ? X *J ^ ijn * P.^ $ x q.J ^ ^ 


t C mnpa = -^ ( o *m.t ] { o X »J ) ( o X p.r ] ( o x a .* 5 Q 2 ) 


Equation (10) can be transformed into a system of simultaneous non- 
linear algebraic equations by division of the volume into an assemblage 
of finite elements and use of isoparametric interpolation 


O 


X h k (°x*) 
A«1 


N * 

- X h k 


(13) 


in each element, where N is the number of nodal points in the isoparametric 
element and h k are appropriate interpolation polynomials. Integration is 
carried out numerically, usually by Gaussian quadrature. Reference [2] 
first linearized Eq. (10) by replacing e e r» by «>•« and 8(©e//l by 6(©*#1 
and then perform modified Newton iterations, setting 


t+ttyM . »+Ar w u-i> + Lu jk) 


in which k is the iteration number and r+At^(°) = t u ^ 



As given in Reference [3] the substitution of Eq. (13) into Eq. (10) yield 
the following matrix equations. 

TABLE 1 FINITE ELEMENT MATRICES 



or 

co 













TABLE 1 (cont'd) 


ANALYSIS TYPE 

INTEGRAL 

MATRIX EVALUATION 

B. MATERIAL 

NONLINEARITY 

ONLY 

f 0 

/ C Ijrs e rs *1) dv 

°V 

* K u = (j B l C b l ° dv ) u 

°V 

jXl 66 11 0dV 

°v 

l F = fi * t °dv 
°V 

C. TOTAL 

LAGRANGIAN 

FORMULATION 

f 0 

J 0^ijr9 O e rs ^O e lJ 
°V 

o k l ° * (fX o c o b l ° dv ) “ 
°v 

fiil ‘o\j ° dv 
°v 

oV “ = (fill, 0 s oV ° dv ) “ 
°v 

ft 0 

Jo S lj fc ®O e iJ dv 

°v 

t_ _ At ti o. 

o F - J o b l o s dv 

°v 






ANALYSIS TYPE 


INTEGRAL 



1 Si i . e 6. e. . Slv 

1 t ijrs t rs t ij 


Sr 

D. UPDATED 


LAGRANGIAN 

|\j 5 Aj * dv 

FORMULATION 

•/ 

Sr 


f* T ij 't e u tdv 


*v 


CT> 

in 
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Numerical Time Integration 


The incremental equilibrium equations must be solved at each time step 
using a numerical integration scheme. In SAP7 the Wilson e - method and 
Newmark method are used. This section follows the formulation given in 
Reference [3]. 

In the Wilson 6 - method a linear variation of acceleration is assumed 
over the time increment t = 6 tt .where (for unconditional stability 
in the analysis of linear systems) e s 1.37, and the equilbrium equations 
are considered at time t + t , 


M 


t + T*. 


+ C t + T U 4 Sc u 


t + T, 


R - 


0) 


where R= R=6( R- t R), and u is the change in displace- 
ment vector during the time interval ttot+r, i.e. u = t+T u - t u. 
Using the linear acceleration assumption it follows that 


t+T. 

U 


t. 

u 


4 


T 

2 


, t4T.. 
( U 


4 



( 2 ) 


t + T t 

U = U 


t. 

u 


2 


t4T.. 

U 


«• . 
U ) 


(3) 


which gives 


t4l„ 


u 



6 

T 


t. 

U - 


0 t.. 

2 u 


(4) 


and 



T 

2 u - 2 


Si 


(5) 


Substituting the relations (4) and (5) into Eq. (1), an equation 
with u as the only unknown is obtained. Solving for u and using the 
linear acceleration assumption, the required displacement, velocity and 
acceleration vectors are obtained. 
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t+At 


• • 


u 


( 1 




t+T« 

U 


( 6 ) 


t+tu 



( % u + t+At ff > 


(7) 


t+ ^u 


t 


u 


+ At 



,t+At.. 
( u 


-t «. . 
+ 2 u ) 


(8) 


EQUILIBRIUM ITERATION 


In order to avoid large integration errors we may choose to iterate 
in each load step until the equilibrium equations are satisfied within a 
given tolorence. For a single finite element, in the TL formulation the 
equation used for iteration is given as, 

< X ' OV > (9) 


i = 1,2,3 ... 


where 


t*At u (l) . t + At u (l-l) + 6 u (i> . 

A more detailed description of the equilibrium iteration equations 
can be found in Reference [3], 
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APPENDIX B 


SAP 7 SUBROUTINE ELPAL 



nnnonnrvian no o non noon noon noon nnno noon noon 
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t •• • •♦•••••* • • ••••• • 

THIS PROGRAM CAN BE USED IN SINGLE PRECISION ON COC COMPUTERS 
AND DOUBLE PRFCISION ON IBM.UN IVAC.DSC VAX 780/11 AND PRI*E 
COMPUTERS. ACTIVATE* DEACTIVATE OR ADJUST THE A30VE CARD FOR 
SINGLE OR 9QU3LE PRECISION ARITHMETIC. 


1ST 
I SR 

SP? 

RATIO 

OELE°S 

DELSIG 


NUM3ER OF STRESS COMPONENTS 
NUMBER OF STRAIN COMPONENTS 
“I r €SSES AT THE END OF THE PREVIOUS UPDATE 
TRAINS AT THE END OF THE PREVIOUS UPDATE 

®ART OF STRAIN INCREMENT TAKEN ELASTICALLY 
INCREMENT IN STRAINS 

I NCREMFNT IN STRESSES* ASSUMING ELASTIC BEHAVIOR 


PROP! 1) YOJNG S MODULUS 
PROP! 2) POISSON S RATIO 

PR0PJ3) INITIAL YIELD. STRESS IN. SIMPLE TENSION 
PROP! 4) STRAIN HARDENING MODULUS 

IPFL = l* MATERIAL ELASTIC 
* ?* MATERIAL PLASTIC 


COMMON /EL/ !ND*ICQL*NT*NPAR(20 )« NUMSG* NSGL *NEGNL * I M A SS ♦ I DAMP * I STA 
1 * NOTE* K5. IN*TFIG*IMASSN.IDAMPN 

COMMON /VAR t NO. KPR I . MO DSX ,K S TE P. ITS* IT=H AX . I REF. I E QR E c * I NOC MD 

^HmS^/VMI SES / Al. 31.C1 .01 ,A?.? 2.C2.02.YL9. 8M.TSR.IST 
COMMON /MATMOD/ STP. ESS ( 4 » , STRA I N 1 4 > , C t 4. 4 ) , THE STR { 4 > .TEMP » I PT. NEL 
COMMON /DtS"^/ DIS0C9J 

J 1 “St I p LT . I PTEL. IPLTDF* 1PLTND ( 2 * 8 ) « 

SNUMEL »NPTR. IPS* I PES 

COMMON /mi.NK/JWR 

DIMENSION PROP(ll«StGfll«EiPSril 

DIMENSION TAU(4).CELSIGI 41 * DELE’ SI 4 ) *05® S I 4) .STATE ( 2 ) 

EQUIVALENCE (NPARI31.INDNL ).(NPAR(5> * ITYP2D » » ( OELEPS *DE°S ) 

DATA NGLAST/ 1000/* STATE /I HE. t HP / 

SQRT(OJ=OSQRT(Q) 

THfr^DGUMV.SVBE’u^SDMrifNGLrS^ClIioN’ofi'eDreDfiPufERS - 

AND DOUBLE PRECISION ON 1 3M* UN IVAC ♦ DEC VAX 780/11 AND PRIME 
COMPlJTS D S. ACTIVATE*D C ACTIVATE OP ADJUST THE ABOVE CARD FOR 
SINGLE OR DOUBLE PRECISION ARITHMETIC. 


IF UPT.NE.l) GO TO 110 
I ST *4 
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c 

8 


c 

c 

105 

C 

110 

c 

c 

c 

c 


TF tITYP2D.E0.2> !ST=3 

IP? ?TYP20. = 0.0) ISR=4 
YM=PR0PU> 

PV=PROP!2> 

Ol=PV/CPV - 1.0001 

«WKtt62WWlfl.ooo-i.ooo*-v. 

C?=PV/ ( 1.000-2 . ODOFPV) 

D2=YKFP»0P<* )/C YM-PR0PC4 ) ) 

§^=Y*/?lJ8Bo - 2.030*PV) /3.000 

jc C ITY , »2D.FQ.2) 00 TO 105 

PLANE STRAIN t AX ISYMMETR IC 
B1=A2*C? 

A1=P1+A2 
GO TO 110 


i 

C 

r 

c 

c 

c 


PLANE STRESS 
A1=YM/(1 .O0-?V*PV> 

31=A1*PV 

YLO. = Y1 ELD 

1. CALCULATE INCREMENTAL STRAINS 

00 120 1=1, ISR 

20 DELFPS C I I = STRAIN! I 1 - EPS € I > 

CHECK FOR WRINKLING BEHAVIOR 
33 333333333 3330303 333 3 33333 303 33 333 33 3333 

J° = l 

C JF( JWR.FQ.OICO TO 30 

EBS= 1 .E-9 

A = ST.RAIN(1)*$TRAIN( ?) 

E=E12*E12 

El=(4+S0PT(?*P+E))/2 
E 2 = C A- SORT C9 s >3+E)l/2 

1 c t El I ?4 » 2C» 23 

24 JP=2 

00 25 1=1,4 
STGC 11=0.30 

25 OELSIGCI)=0.000 
GO TO 150 

20 EP=PV*E1+E2 

IFCEP+EBSI2&*30,30 

26 JJ>=3 
P=3/CF1-E2> 

Q=E12/tFl-E2) 

Pl = ( 1+P >*Y*/2.00 
P2=( l-PIFYM/2.00 

Q1=3*YM/4.D3 
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8 

c 

c 

i 

30 


s'? = YM/ 4 

6?LSIGCl)*( , l*0FLE°Sfll , * , 01*0FLE o SI3) 

S!G< l)sP1*Ei>S(l) + Ql*E?SC3 » 
DElSIGi2>=P?*OELe?Sm«-91*DEL5PSm 
SIG<2lsP2P€PSC2l+01PEPSC3> 
DEL^!G(3l*0tPC0FLEPS(l) + 05L.FPS{2>)+02*0ELEPSC3l 

0Eafmmrs5G sn,+EPS(2M " Q? * EpsC3> 

SIG(4)=0.*>D0 
GO TO 150 

33 3353333333333 3339333 33333333333333333333933 


C 

C 

C 


150 

£60 


C 

c 

c 


c 

c 

c 

c 

c 

170 

150 


C 

c 


CALCULATE 

ASSUMING 


THE STRESS INCREMENT, 
ELASTIC BEHAVIOR 


♦ 

♦ 


91PDFLEPS (?) 
A1*DELEPS!2I 


D e LSIG(n * AlPOELEPSm 
OFLSjsm » 3i*0FLEPsm 
D5LSIGC3) * C1PDELEPS(3> 

O e LS I G ( 4 1 a 0,300 

IF ( ITYP20. c 0.2) GO TP 150 

DELSISIM =_81 P I DELE PS 1 1 > ♦DELE PS 12)) 

I e ( ITYP20, = '3.1 ) GO TO 150 

0FLSIGC1) = DELS IG( 1 ) ♦ 31PDELEPSI4 ) 

DELS TGI?) = OFLSISm ♦ B1 PDFL EPS! 4 I 
DELS I G ! 4 ) = OELSIGC 4) + A1PDEL EPS! 4) 

3 « CALCULATE TOTAL STRESS C S, 

ASSUMING ELASTIC 8EHAVI OR 

TAUI4) a 0.000 
00 160 1=1, 1ST 

TAum = siscn + oelsigid 

4. CHECK WHETHER PTAUP STATF CF STRESS FALLS 
OUT SI OE THE LOADING SUR = ACE 

SM = !TAU!l)+TAU(2)+TAUf 4)1/3. 

SX a TAUI1) - SM 

SY = T4JI2) - SM 

SS = TAum 

S7 = TA’J ( 4 ) - SM 

IF( JWR.EQ.11G0 TO 170 

FT =.SQ0P(SXPSX+SYPSY+S7*SZ) ♦ SSPSS - YLDPYLD/3.D0 
IF (cj| 170*170,300 

STATE O e STRESS WITHIN LOADING SURFACE - ELASTIC BEHAVIOR 
IPEL=t 

ST^E S S ( 4 ) a 0.000 
DO 130 I = 1 , I S T 
STRESS ( I ) = TAUCI) 

IF ! ITY' J 2D. p Q.2 ) STRAIN! 4) = EPS ( 4 I ♦ D!P( DELEPS ( 1 ) ♦ DELEPS(2)> 
T F f J°.GF.2)STRAIN!4)=EPSI4) 

GO TO 400 
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STATE OF STRESS OUTSIDE LOADING SURFACE - rL ASTI C BEnAVi OR 
IF (IPEL.EO.l) GO TO 320 


iWAS ®L ASTIC 


4SfV5 2 - o. 

DO ?15 1=1 ,T ST 

tauu) = sism 

GO TO 370 


.WAS ELASTIC 


DETERMINE PART OF STRAIN TAKEN FLASTICLY 

i I °EL = 2 
IPEL=2 

SM = (STG(1)*SIG(2) + SIG(4))/3.0D0 
SX = sic.m - SM 
SY = SIGm - SM 

i? : IfBJI! - s* 

DM = (DELSI3(l)+DFLSIS(2 )+DELSIG(*))/3.0D0 
OX = DELSIS(l) - DM 
OY = OFLSIGC 2 1 “ DM 
DS = DFL SI 5C 3 ) 

DZ = DELSI3I <► I - OM 

A = OX*DX OYSPY 2.0D0*DS*TS + DZ*OZ 
3 = SXSDX ♦ SY^DY ♦ 2.0D0^SS*TS ♦ SZ*OZ 

E = SX*SX ♦ SY*SY + 2.0D0*SS=SS ♦ SZ*SI - ?. . OCO*YLD=YLD /3 . ODO 
RATI0=(-B*S?RT(3*3-A*E)) /A 
O'* 350 1=1, 1ST 

TAU (!) = S 1 3 ( I ) ♦ RATI3*0ELSIG( I) 

IF ( ITYP2D.EQ.2 ) STRAIN! 4)=EPS (41 ♦ R ATI0«D1*( DELEPS (1) 

1+ D C LEP3 121) 

*TAIJ* NOW CONTAINS (PPEVIOUS STRESSES ♦ 

STRESSES DUE TO ELASTIC STRAIN INCREMENTS) 

5. CALCULATE PLASTIC STRESSES 

DFTERMINE INCREMENT INTERVAL 
M=2O.DD*S0RT(eT)/YL0*l 
1* (M.GT.30) m=30 
XM = (1,00 - RATIO ) /M 


DO 380 1=1 , 1 SR 
DE°S(I1 * XNFDELEPSU) 


CALCJLATION OF EL ASTOPLASTIC STRESSES 


(START) 
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4 


d?ooo 

c 

c 

c 

C?001 

C 


? 6 ° 

C 

c 

c 


r 

£ 

c 


c 

c 


c 

590 


C 

£00 

c 

c 

3°0 


IM5ttfci??SifPS! 


rfO.ITE(6*FI 


I DEPSI 11,1=1,4) 

WR!T=f6»*l A,B,S, RATIO 

WHITE 16*2091 ) 

FOP. MAT f ?X, 3HT AU) 

WRITE(6»F) ITAUd),I*l,4l 

00 600 !.W=1,M 

CALL MIDEP f TAU,DEPS,C) 


DO 560 1=1,1 ST 
03 560 J=1,ISR 
TAUdl = TAJ(!» 


CORRECTION 


♦ Cd,J) * OEPS(J) 


2" * fTAUd)*-TAU(2)+TAUt 411 /3. DO 

OX = TAUdl - DM 

OY = TAU(2I - 01 

CS = TAUC3I 

07 = TA’J(4I - OM 


IF (PR0P(4).E?.0.) 00 TO 530 

strain-hapd-ning MATERIAL 


UPDATE YLD 


0915=1.503 

YL0=S9RT f 3P159(DX9DX*DY*DY+2.*DS*DS+0Z*0ZI 1 
^0 600 


GO 


PERFECTLY PlASTIC MATERIAL 
FTA=.509*(DXRDX + DY*DY + DZ*3Z) ♦ DS*DS 


FT= C T A - FT3 


TO 600 


IF I FT.FQ.O) 53 ... ... 

IF ( ITYPZ0.EQ.21 GO TO 590 

Core-. i . do i-SCRT(FT9/FTA ) 


TAUdl = TAJdl 
TAUI2I = TAJC2) 
TAU(31 = TAJC3) 
T AUt 4 1 =TAU( ^ 1 + 
GO TO 600 


C0EF*0X 
C0EF$DY 
♦ C0EF*QS 
C0EF9DZ 


C0E c = SQRT( FT3/FTA) 
TAU( 1 ) =TA'Jt l l^COEF 
TAU(2)=TAUdlPCQEF 
T AiJ( 3 1 =TAUt 3 1 EC0EF 
STRAINt41=STRA!N(4) 

CONTINUE 


♦ (COEF - 1.00 )*0M/BM 


CALC JLATION OF EL ASTOPl ASTI C STRESSES 


STRESS (4) = O.DO 
00 390 1=1, 1ST 
STRESSCI) = TAUdl 


t END 1 
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£ 

c 

400 

410 

420 


6. UPDATING STRESSES* STRAINS* YIELD. NS 

00 410 1=1. 1ST 
S I G C I > = STRESS!!! 

DO 4 JO 1=1, ISR 
E°sm = STRAINCM 
YIELD = YLD 

IF C ITY°20. = 0.2> EP S ! 4 ) = STR AI N ( 4 ) 

IF (KPRI.E3.0) GO TO 700 


£ 

C 

c 

c 



4*0 

460 

£ 

C 
C 

c 
c 


IF ! ICCUNT.EQ.3) RETURN 

7 # . FORM THF MATERIAL LAW 

IF { I^EL.FO.l) GO TO 450 

ELASTO-oLASTIC 

CALL MIOEP I TAU*DEPS,C ) 

WRITEI6.300?) 

FftRMAT(5X.l5HTAU,D c PS,CI I.JI) 
WRITEI6.P) “ 

W* I r E! 6* P ) 

WRITE! 6.P ) 

RETURN 


f TAU(II*I=1 *4) 
(0EPSm,I=l,4) 

( ICII, J),J=1*4)*I=1,4) 


ELASTIC 

00 460 1=1, ISR 
00 460 J = 1 * l S R 
C! I, Jl=0.00 


MODIFICATION 
FOP. WRINKLED 


OF STRESS-STRAIN MATRIX 
BEHAVIOR 


33333333333333332339333 33333333333333 
IFCJP.NF.2) GO TO 36 
rn 40 i=i,3 
O'T 40 J=l,3 

40 C( I , J)=0.000100 

R C TURN 

36 I C (JP.NE.3)*;0 TO 35 

c(i,n=pi 

C 1 1 * 2 ) =0,003 
Cl 2,1 >=0.020 
C(l* 3>=01 
C ( 3, l ) =Q1 
C( ;, 2)=P2 
CC2, 31=01 
C I 3, 2 I =01 
C ! 3, 3 ) =02 
RETURN 

C 333 3 3333 3333333 33393333 33333 3233 

35 C(1*1)=A1 

C! 2,l)= a l 
CCl.2)=31 
CC 2* 2 )=A1 


4 
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IIypIo.NE.O) RETURN 

ip uw&muvwiio 


C(1.4)*Bl 
C(2.4)*31 
C(4,l)=81 
C ( 4. 2 ) =B1 
C(4,4 )*A1 

RETURN 

C ( 4. 1 ) =3 2 
C 1 4. 2 ) -3 2 
C(4.3)=0.C 
C(4,4)=A2 

RETURN 


P R I N T I N 


OF STRESSES 


IF CIPEL.E9.1) GO TO 705 

PM=( STRESS! 1 1 ♦ STRESS ( 2 ) + STRESS ( 4) ) /3 .0D0 

DV=STP.ESS(1) - OM 

0Y= STP ESS ( 2 ) - OH 

DS= STR £S Si 3 i 

P7=STPESS(4I - DM 

FT=.5D0*(0X*DX ♦ DY*DY + 02*02) «■ DS*DS - YLD*YLC/3. 000 
I c (INDNL.NE.2l GO TO 500 

IN TOTAL LAORANGIAN FORMULATION. 

CAUCHY STRESSES ARE CALCULATE") AND PRINTED 

CALL CAUCHY 

IF (NG.NE.NSLAST) GO TO 802 
IF ( NEL.GT.NELAST) GO TO 506 
IF (IPT-1) eiO, 803.810 

NGLAST = NG 

IF ( ITYP2D) 303,305.803 

.jwrfiyisjii 

GO TO 306 
IF( IPS.NE.O) 

*WRI T E C6, 20031 

nflast=nel 

Ir(IPS.NE.O) 

3WRITE (6.2004) NEL 
CALL MAXMIN I STRESS .SX.SY. SM) 

IF ( ITYP2DI 313.815.813 

IF ( JWR.EO.l ) GO TO 817 
IF (IPS.NE.O) 

4WRITF (6.2005) l PT, STATE ( IP FL) , 

1( STRESS! I). I =1.3), SX.SY. SM. FT 
GO TO 818 
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317 

820 


813 


C 

315 


8151 

C 

C 

200 ? 


20 03 


20 C4 
20C5 

mi 

2003 
20 0 7 
C 

2010 


J (6f25lS) N ?£^itAf§?}pEl )«( STRESS!! )» 1*1.3) *$X*SY*SM 
*JoiTE(I^200?l 4 lPT*C STRilstI>*I = l*3)*SX*SY* SM 

i*'?f p??.ne.o) 

PtfRITEC 13*2210 ) TIME.NG.NEl.IPT.I STRESSCI )*I=1*3).SX.SY* 
1SM.«=T*JP 
RETURN 

*W?1TE S ( 5^2007) !PT» STATE ( IPFL) * STRESSC4) * 

1C STRESSC I).I = l*3).SX.SY. SM* FT 
IP ( IPES.nI.O) 

*?WR ITEC 18*20101 TIME*NG*NEL* 1° T. STRESS! 4 )*( STRESS C ! )• 1 = 1*3)* 
23X.SY. SM.FT, J° 

CONT I NUE 
RETURN 


STRESS 

max Stress 


L :MENT STRESS 
STP.ESS-YZ 

33H NUH/IPT STATE 


FORMAT ( 50M ELEMENT 
14PMSTRESS-YZ 
25 n H NMM/I»T STATE 
34<M 

4EHE*JNCTIQN / 

FOR mat (* 0-f P 

1 54P STRESS-? ? 

221X.5HYIELD' 

345H 
43-H 

FORMA T f It/) 

FORMAT f 5X * T 
FOR ’RAT (EX* 

FORMAT ( 5 X * I _ _ 

c ?RMA t C5X.I 2*2X,7HWR!NKIE. 1X.3F 
FORMAT ( 5X * I “ " *“ 


STRESS-YY 
MIN 


STRESS-ZZ 
TRESS»21X*5HYI ELD/ 


ANGLE. OX* 


STRESS-XX 
MAX STRESS 


STRESS-YY 
MIN STRESS. 

t 

ANGLE. 9X* SHF'JNcflON/) 


!2*2X*Al* 6HL ASTI C* 1 X. 

2*2X*7HWP. !NKLEtlX*3F14«6*3X*2 f: 14«6*3X*F6»2) . 

2 » 2X, A 1 *6HL A STI C * 1 X* 4G 14.6.3X* 2S14.4 *3X*F6.2 *3X* El 4 


3E 14.6 *3X. 2E14,6*?X*F6.2 *3X. El 4 

*3E14.6*3X.2E14.6.3X*=6*2> 

14«6*3X*2E14«6*3X*F6»2) 


FORMATfFin.5 ,315* 5**3E14.5/ 5E14. 6.2X.I2) 
END 
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APPENDIX C 

SAP 7 INPUT LISTING FOR VERIFICATION EXAMPLE NO. 1 


-I 
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***Tf;ST P9Q3LEM FOR 
»CCNTqcL 
Ml *6Q,NL»E 
M4, 0,0, 321 
*5,320,0.0625,0.0,4 
*L0£0S 
3,2,3 


NASA PROJECT 


f?5 


0.0, 0.0, 1.0,0. 04, 20.0,0.04 

2,3 

0.0,0.0,1.0,0.0,20.0,0.04 
L, 192,192, 10,2,1,0.003 
L, 193,191, 1,2, 1,0. 006 

1.11.131.170.3.1.0. 002 
L, 17, 170, 17, 3, 1,0.003 
L, 29, 164, 17, 3, 1,0.004 
L, 1,1 71, 170, 3, 1, -0.002 
L, 12, 165,17, 3, 1,-0.009 
L, 18, 154,17,3, 1,-0.004 
L, 192, 19 2, 1,2, 2, -0.014 
L, 13 3, 193, 1,2,2,-0.024* 

L, 134, 1*4, 1,2,2, -0.019 
L, 195,135,1,2,2,-0.012 
L, 196, 186, 1,2,2,-0.006 
L, 193, 199, 1,2, 2, 0.006 
l,199»l q ?»l, 2*2,0. 012 
L, 190, 1*0*1,2,2,0.013 
L, 191, l«>i,l, 2,2, 0.024 

I, 192,1*2,1,2,2,0.014 
*2/0, NL 

J, 2, 3, 2, 0,9,1 

7. 2E + 4, 0.333,1.09*4 
*9FA M» L 

9.0. 1 ,1,8 

9.0E+6,0,3.0E+6,0» 0 

1.0. 0.0.0.0*0.1666667,0.0333333,0.0333333,0.166667,0.166667 

SIItch 

FINISH 
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*«TF$T PROBLEM FCR NASA PROJECT! MODEL INPUT 

192.71 

ECHO 

CC START GENERATION 0*= NODES 

NOCTtltO. 0*0.0 *0.0 

TS *10*1*1*1*1* 0.0*0.0*0.6 

T9,l 0.1* 11 *1*17*0. 0*1. 2 *0.0 

NDOS* 1 2* 0. 0 *0. 5*0. 0 

TR, e , 12* 12* 1*1*0. 0*0. 0*1. 2 

TR, 7,12*17, 1,1 7, 0.0 *1.2 *0.0 

N?DO. lSf.0.0.1 2.1 *0.0 

T'** 10* 152*182* 0*1* 0.0. 0.0*0. 6 

CC START GENERATION OF 3CUNDAR Y CONDITIONS 

FIX. 1.1.1*02 

?!X»1*11*1*DY 

FIX, 1*181*1* DX» NR 

FIX, 192*192*1* DX,RY*RZ 

CC SAVE NODAL CO-ORDINATES INFORMATIONS 

cc r ^Tart GENERATION OF ZD ELEMENTS 

FL *1*20*1*0. 0*0.01*1*8 

l,l r . *?0*3*12*15*13*2 

PEP. d, 1 , 1 , 1,17 

FL.l 1. 2D, 1,0. 0,0. 01, 1,8 

3* 20*22*5, 13, 21, 14, A 

REP, 9, 11, 11,1*17 

FL.21 ,20.1,0.0,0.01,1,8 

5,22, 24,7,14,23. 15,6 

REP, 9, 21*21* 1*17 

EL, ■*1*20,1, 0.0*0. 21, 1,3 

FL, 41, 20, 1,0.0, 0.01, 1,8 

9,26,23,11,16,27,17,10 

R , rP,9,4t,*-l *1,17 

CC SAVE ELEMENT INFORMATIONS 

EWRITE 

CC START GE.NE’ATUN OF BEAM ELEMENTS 
FL,1, BEAM, 182, 183* 103,1,0,0,0,0.0,1.00 
EL,2, B, 133, 134,104,1,0,0,0,0.0 ,1.00 
FL ,3, 5,134,135, 105, 1,0, 0,0 ,0.0,1. 00 
EL ,4 *5*195,1 36 ,105*1,0*0,0, 0.9, 1 # 00 
EL, 6, 3,186, 1 87, 107, 1, 0,0, 0,0.0 ,1.00 
EL ,6,3 ,197, 13? ,103, 1,0, 0,0, 0.0, 1.00 
EL, 7 ,3 ,193, 189, 10? *1,0, 0,0,0.0,1.00 
EL ,3,5,139, 190,110, 1,0, 0,0, 0.0 ,1.00 
FL,9, 3, 190, 101 ,m ,1,0,0, 0,0.0 ,1.03 
FL,1 0,3, 191 ,192,11 2, 1*0,0,0,0.0,1.00 
FVJR1 T? 

CC S T A°T GENERATION OF CONSTRAINTS FLEM^N^S 

FL,1 »CP,1 

171,182 

REP, 4, 1,1, 1,1 

EL, 6, CP, 1 

177,188 

REP, 4, 6, 6, 1,1 

ENRITE 

FL,1,CP,1 

176* 1 87 

EWRITE 

END 


INFORMATIONS 
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APPENDIX D 


SAP 7 INPUT LISTING FOR VERIFICATION EXAMPLE NO. 2 
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MEHe,?ANS( CIkCUMF - ,AnC * a >- TP ANu SPJFDR NASA PROJECT 

M 1 *20 * NL * F 
M4, 9*0* 241 
K5»H0t 0. 0625*0*0*4 
*LCA0S 

0*1*3*0*0*0*48 

1*3 

0. 0*0*0 *1.0 *0.02*1 5.0 *0.02 

2*3 

0*0*0.p*1.0*0*0»15.0*0.02 
P29*7, 19,1 1,0* 2*1* 2*0.0*0.5*0.01 
P20* 19, 29, 2 2,0,2*1,2*0.5,0.666025*0.01 
P20* 29, 40, 3 3*0, 2. 1.2* 0.66 6025* 1.0,0.01 
P‘ , 0*4p,5l, 4^ *0*2*1*2* 1.0,0.966025,0.01 
P20. 5 1*62 *55, 0,2* l *2, 0.866025*0.5*0.01 
P?0*62*73*66*0*2*l *2* 0.5 *0.0«0.01 
P2D, 73, 94, 77, 0,7, 1,7, 0.0 *-0.5. 0.0 l 
£20,9 4 *05,89*0.2*1 *2, -0.5.-0.8660 25 *0.01 
P?9»95*1CA»98*0*2* I * 2 * — 0. 966025*— 1.0*0.01 
P7D* 106.117* 110*0, 2* 1 *2*- 1.0,-0.966025 *0.01 
P?0,117*123*121*0» 2,1 , 2.-0. 866025 *-0.5 *0.01 
P’D* 12 9*7, 13 2, 0,2, 1,2, -0.5, 0.0. 0.01 
P?D,7 ,18*1 1,0, 1*1, 2, 1.0*0.866025* 0.01 
P20* l 8, 29, 2 2, 0,1 ,1,2* 0.96 602 5* 0.5 *0.01 
P23, 2°, 40, 23, 0,1,1, 2, 0.5 *0.0,0.01 
£*£»£? *?♦ 0.0 *-0.5* 0.01 
P 20*51, 62*55,0.1,1 ,2, -0.5 ,-0.966025* 0. 01 
P20*62, 7 3,66,0,l,l,2*-0.866025 ,-1.0*0.01 
P7D, 73, 4 4, 77, 0,1,1 ,2,-1. 0 *—0.3 660 25,0. Cl 

P20, 95, 105, 09, 0,1,1, 2,-0. 5*0.0 *0.01 
P?0* 106 *11 7, 11 0,0, 1,1 *2, 0.0 *0.5,0. 01 
P 2D, 11 7, 12 9, l?l,0, 1, l ,2. 0.5,0. 86 6025,0.01 
J*2D»128. 7 .132*0*1*1, 2*0.56 60 25, t.O, 0.01 
*»’C,7, 18*1 1*0* 2* 7* 2 ♦■‘1.0, -0.96 6025* 0.01 
P 20, 18* 29* 2 2 *0,2, 2 » ?*-0. 8 66025,-0.5*0.01 
P?D, 2 «» *40, 3 3, 0,2*?, 7 *-0.5 ,0.0*0.01 
P20.4A* 51, 44, 0,2, 7, 2 *0.0, 0.5, 0.0 l 
P^C, -si ,67, ,9,'>,2,2* 0.5, 0.86 6025 *0.01 

P7D* 5 2,^3 *6 6, 0,2* 2* 2*0.96 602 5, 1.0 ,0.01 
P2D, ■'3, ?4, 77*0*7, 2* 7, 1.0, 0.356025* 0.01 
P?0» 9 ^*f5*9®,0* 7, 2, ?, 0.96 6025* 0.5,0.01 
P 20*95,1 06,99, 0,2, 2*2* 0.5 *0.0* 0.01 
P?n,iot,ii7, 11 0,0, 2, 2 *2. 0.0 *-0.5, 0.01 
?20, 117, 129 , 121,0, 2*2 ,2»-0.5,- 0.366025*0.01 
P? D* 1 2 8 » 7,132*0,2*7,2 *-0. 86 602 5* - 1 . 0 ,0.0 1 
°2D, 7. 19, 11, 0,1,2,2,0.0*0.5*0.01 
P’O, l «, 79,22, 0,1, ?, 2, 0.5, 0.866025 *0.01 
°?0» 2° *40, 33,0,1, 2, 2,0.26 60 2 5, 1.0,0.01 
P2D, 40, 51, 44,0,1,2, 2, 1.0,0.966025*0. 01 
P 2D, ;:1 ,62, 5 5,0,1 ,?, 2* 0.866025, 0.5 ,0.01 
P.7P, 62, 73, 66, 0,1, 2, 2, 0.5, 0.0, 0.01 
P7D* 7 ?, 34, 77, 0,1, 2,2,0. 0,-0. 5, 0.01 
P2 0,34,95, 8 p ,0,1 ,2 ,2,-0. 5 ,-0.966025,0. 01 

P2D,o5,l06,og,0,l,2,2*-0.966025*-1.0,0.01 

P 2D, 106 *117, 11 0,0, 1,2* 2, -1.0, '-0.866325,0.01 
P2D, 117 ,123, 12 1,0, 1 ,2 ♦ 2 » — 0 • 966025 ,-0.5 *0.01 
P|9,12 3*7*132,0*1, 2* 2 *-0.5,0. 0,0. 01 

1 »2.S*2*C»9»1 


7.7=+4,0.333,1.0E44 

FINISH 
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CSC! OCULAR ME^BRA'JFCCIRCUMFERANCI AL 

1*2*36 

ECHO 

CC GENERATE NCO c S 
N0DE*l,0.-),20.0.0.0O 
TR,6,l,l*l* 1*0.0*10.0*0*0 
RO T * ltl.lt 1,7*-15*0,X 
R“)T»l,3,3*l*6,-i5.0,X 
RQT,i,5*5*l,5*-15.0*X 
R?T, 1*7,7, 1*4, -15. 9.X 

CC GENERATE 8. CONDITIONS 
ETX.l ,122,ll,DY»ni 
FIX, 3*129*11, DY. 01 
V.VP I Tr 

CC GENERATE ELEMENTS 
p LE*l, 20*1*0*1. OE-2*l*8 
14,12*1,3* 13*8*2*3 
EL5. 2*20*1, O.l.OE-2, 1,3 
16,14,3,5*15,9,4,10 
ELE, 3,20, 1*0,1 *05-2* 1,8 
1?, 16,5,7,17,10,6. 11 
RE, 10, 1,3*1*11 
EL.34, 20, 1,0,1.0E-2*1 ,R 
3, 1, 122,124, 2* 129*123, 130 
EL*3S,2D*l,0,1.0r-2*l,3 
5,2,124*126*4, 170,125,131 
FL,7S* 20,1 *0,1. HE- 2, 1,8 
7, 5,126*123*6* 131*127,132 
ENRITE 
END 


AND S) FOR NASA 




t 





1. Report No. 

NASA CR- 172325 


2. Government Accession No. 


3. Recipient's Catalog No. 


4. Title and Subtitle 

Finite Element Analysis of Wrinkling Membranes 


5. Report Date 

April 1984 


6. Performing Organization Code 


7. Author(s) . 

Richard K. Miller, John M. Hedgepeth, Victor I. Weingarten, 
Prasanta Das, and Shahrzad Kahyai 


8. Performing Organization Report No. 

USC-CE-83-05 


10. Work Unit No. 


9, Performing Organization Name and Address 

Dept, of Civil Engineering 
U. of Southern California 
Los Angeles, CA 90089 


11. Contract or Grant No. 

NAGl-235 


12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, DC 20546 


13. Type of Report and Period Covered 

Contractor Report 


14. Sponsoring Agency Code 


15. Supplementary Notes 

Langley Technical Monitor: Dr. Wilbur B. Fichter 

Consultant, Astro Research Corporation, Carpinteria, California 


16. Abstract 


The development of a nonlinear numerical algorithm for the analysis of stresses 
and displacements in partly wrinkled flat membranes, and its implementation on 
the SAP VII finite-element code are described. A comparison of numerical results 
with exact solutions of two benchmark problems reveals excellent agreement, with 
good convergence of the required iterative procedure. Also reported is an exact 
solution of a problem involving axi symmetric deformations of a partly wrinkled 
shallow curved membrane. 


17. Key Words (Suggested by Author(s)) 


18. Distribution Statement 


Membranes 
Wrinkl ing 
Finite element 


Unclassified - Unlimited 




Subject Category 39 

15 Security Classif. (of this reportl 

20. Security Classif. (of this page) 

21. No. of Pages 

22. Price 

Unclassified 

Unclassified 


85 

A0 5 


For sale by the National Technical Information Service, Springfield. Virginia 22161 




















t 




